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Dissipative processes in non-equilibrium many-body systems are fundamentally different than 
their equilibrium counterparts. Such processes are of great importance for the understanding of 
relaxation in single molecule devices. As a detailed case study, we investigate here a generic spin- 
fermion model, where a two-level system couples to two metallic leads with different chemical po- 
tentials. We present results for the spin relaxation rate in the nonadiabatic limit for an arbitrary 
coupling to the leads, using both analytical and exact numerical methods. The non-equilibrium 
dynamics is reflected by an exponential relaxation at long times and via complex phase shifts, lead- 
ing in some cases to an "anti-orthogonality" effect. In the limit of strong system- lead coupling at 
zero temperature we demonstrate the onset of a Marcus-like Gaussian decay with voltage difference 
activation. This is analogous to the equilibrium spin-boson model, where at strong coupling and 
high temperatures the spin excitation rate manifests temperature activated Gaussian behavior. We 
find that there is no simple linear relationship between the role of the temperature in the bosonic 
system and a voltage drop in a non-equilibrium electronic case. The two models also differ by 
the orthogonality-catastrophe factor existing in a fermionic system, which modifies the resulting 
lineshapes. Implications for current characteristics are discussed. We demonstrate the violation 
of pair-wise Coulomb gas behavior for strong coupling to the leads. The results presented in this 
paper form the basis of an exact, non-perturbative description of steady-state quantum dissipative 
systems. 

PACS numbers: 03.65.Yz, 05.60.Gg, 72.10.Fk, 73.63.-b 
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I. INTRODUCTION 

Over the past several decades tremendous effort has 
been put forth to understand the dynamics of a small 
quantum entity coupled to a thermal bath [l[ . Important 
problems that can be distilled to this form include the 
interaction between localized magnetic impurities and 
itinerant electrons (the Kondo problem) [E Hj] , electron 
transfer in aqueous environments and proton tun- 
neling in biomolecules [E [f| . The study of such quantum 
dissipative systems cuts across traditional disciplines and 
impacts fields from biology to quantum information the- 
ory [H. 

Our detailed understanding of quantum dissipative 
systems is essentially confined to problems which involve 
a single thermal reservoir [7|. In this case, traditional 
measures of dynamical interest are equilibrium correla- 
tion functions or simple measures of the decay of one-time 
quantities when the initial condition is not one of ther- 
mal equilibrium for the global system. While important 
for the understanding of various experimental situations, 
this latter form of non-equilibrium behavior is well un- 
derstood and generically takes the form of an asymptotic 
exponential decay to the thermal equilibrium state or the 
ground state (at zero temperature) [j], 



A less well understood type of non-equilibrium behav- 
ior may manifest when a small quantum system is cou- 
pled to more than one reservoir [E EE EE EE EE EH 
EE EI El El, El, HE US HI- Here, the generic situa- 
tion is one of a non-equilibrium steady state, regardless 
of the initial preparation. Given the fact that this multi- 
bath scenario is standard for prospective single-molecule 
devices [IE [Hj] as well as more general problems, it 
is imperative to understand the fundamental relaxation 
motifs that emerge in such nontrivial non-equilibrium 
cases. Recent work raises the question of whether stan- 
dard tools borrowed from typical equilibrium quantum 
dissipative systems are useful in the steady state non- 
equilibrium case. For example, the simple equivalence 
between bosonic and fermionic baths (as obtained via 
bosonization [IE HE H3|) is lost in the multi-bath case, 
while mean-field approaches are fraught with danger due 
to the fact that a voltage bias may assist tunneling even 
at zero-temperature, rendering the meaning and stability 
of Hartree-Fock minima unclear [10, 123, 129L SO, Ell • 
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While several recent papers have taken up the task 
of describing the steady state, non-equilibrium dynam- 
ics in different model problems, our goal here is a first 
step towards a detailed and systematic understanding 
of dissipative relaxation in the simplest model prob- 
lems resulting from coupling a small quantum system 
to several baths, namely generalized spin-boson mod- 
els [E EE Eli • It should be noted that the term "spin- 
boson" is a misnomer; the interesting and relevant case 
is that of fermionic reservoirs, which dramatically differ 
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from the case of bosonic reservoirs when the system in- 
teracts with more than one bath with different chemical 
potentials. On the other hand, as in the standard spin- 
boson model, it is the physics of the x-ray edge singularity 
[I2> HH HH HH that forms the fundamental building block 
of the description of dynamic observables. Here, it is 
the recently studied non- equilibrium x-ray edge problem 
HE S3, H, SI S3, 5H, SI that lies at the core of the re- 
laxation behavior of standard correlation functions. The 
more complex physics of the non-equilibrium edge behav- 
ior allows for a richer range of dynamical behavior than 
in the well-studied equilibrium case. 

In this paper we will confine our discussion to calcula- 
tions that are perturbative in the bare tunneling matrix 
element of the system, but allow for arbitrarily strong 
coupling to the leads. We will employ both analytical 
and numerical techniques to describe the dynamics. The 
numerical approach involves a computational solution of 
the non-equilibrium x-ray edge problem that is numeri- 
cally exact on all relevant time scales. This will allow us 
to describe the full cross-over behavior from the regime 
where equilibrium effects dominate, to that where the 
full non-equilibrium behavior (such as bias-induced de- 
phasing and complex phase shifts) is manifested. This 
is crucial, since the full frequency dependence of relax- 
ation rates and generalized fluctuation-dissipation ratios 
depend on the entire time history of the dynamics (llj . 

We will demonstrate that interesting behavior oc- 
curs in specific parameter regimes that lead to anti- 
orthogonality effects and bias-induced tunneling. In 
particular, the bias-induced tunneling regime at zero- 
temperature may display a very broad Gaussian decay of 
the polarization at strong system-leads coupling. In this 
regime, the relaxation behavior shows interesting simi- 
larities to the usual high-temperature Marcus (or semi- 
classical polaron) behavior [43, 0, with potential 
bias playing the role of temperature, although crucial 
differences exist that make these analogies imprecise. 
Lastly, we investigate the crucial question of the accu- 
racy of the pair- wise Coulomb gas decomposition for non- 
equilibrium steady state systems. We note that the meth- 
ods discussed in this work form the basis of a numerically 
exact path-integral description of quantum dissipation in 
such non-equilibrium problems [461 ]. 

This paper is organized as follows. In section II we 
describe our model system (the out of equilibrium spin- 
fermion model). Section III presents an overview of 
the analytical results for the non-equilibrium dynamics, 
along with the relation to the non-equilibrium x-ray edge 
problem, while section IV presents numerical results. In 
section V we present the implications for the tunneling 
rate. Our results imply a breakdown of the Coulomb gas 
picture at intermediate times, described in section VI. In 
section VII we conclude. 



II. MODEL 

Our model system consists of a biased two state sys- 
tem (spin) coupled to two electronic reservoirs n = L, R 
held at different chemical potentials. In what follows 
we assume that the temperature is zero, and investigate 
the possibility of voltage activated excitation between the 
spin states. The extension to non-zero temperature is 
straightforward, both analytically and numerically. The 
total Hamiltonian is the sum of three terms: 



H = H S + H ( B f) 
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(1) 



The spin system H$ consists of a two level system (TLS) 
(creation operators d±) with a bare tunneling amplitude 

A and a level splitting B. The reservoir term Hg in- 
cludes two non-interacting metallic leads n = L, R, where 
a non-equilibrium state occurs when the leads have differ- 
ent chemical potentials Ap — /i^ — [ir =/= 0. The system- 
bath interaction H^g couples the spin with scattering 
processes inside the leads (diagonal coupling) , and in be- 
tween each lead (nondiagonal coupling), and we choose 
conventions such that only one of the spin levels couples 
to the leads, 



H s = 



H 
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Here rc^i) = (7 ± cr z )/2 is the number operator, with 
I as the identity operator (47[. The operator a\ n (a fe .„) 
creates (annihilates) an electron with momentum k in the 
rt-th lead. In this paper we focus on the model presented 
in Ref. [36j, where the momentum dependence of the 
scattering potential is neglected. System-bath scattering 
potentials are then given by V n ^ n ' , where n, n' — L, R are 
the Fermi sea indices. Our main conclusions, however, 
are valid for more general cases. 

We assume that the reservoirs have the same density of 
states p(e), typically modeled using a Lorentzian function 



Pie) 



D L /2 
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where Dl is a bandwidth parameter. We typically work 
in the limit of wide bands, Dl 3> Afj,, therefore to a good 
approximation p(/ii) — p(pr)- 

Note that we ignore the spin degree of freedom of the 
reservoir electrons in our discussion. In what follows we 
refer to the energy difference B as a magnetic field, in 
order to distinguish it from the voltage bias A/x. We also 
define two auxiliary Hamiltonians H± that will be useful 
below 

H ± = ±^ + H^[n d = ±]+H { J\ (4) 
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Explicitly, H± includes the electronic reservoirs and 
system-bath interaction, given that the subsystem is in 
the ± state. The model ((2]) contains much of the physics 
of the Kondo model Q, while lacking direct coupling of 
the reservoir degrees of freedom to spin-flip processes. It 
also contains the spin-resonant-level model of Ref. [ll| 
with a particular choice of system-bath couplings. We 
discuss the spin-resonant-level model in more detail in 
Appendix A. 

Crucial parameters of the model are the L and R scat- 
tering phase shifts. In equilibrium the phase shifts are 
given by [H, [3(| 

tan5± = i[(ai + a 2 ) ± [(a x - a 2 ) 2 + 4|i/| 2 ] 1/2 ], (5) 



An important quantity that will be useful below is the 
sum of the phase shifts squared, 7 = — (<5 2 + Sj^/ir 2 . 
While for the general model Eq. ([7|) yields complex num- 
bers, when the system is symmetric (a± — a 2 ), the phase 
shifts are complex conjugates and 7 is real, although pos- 
sibly negative. We discuss the implications of this result 
in section IV. C . 

For the sake of completeness and comparison the equi- 
librium spin-boson model is discussed in Appendix B. 
In the nonadiabatic limit, this model yields the classical 
Marcus rate at high temperatures when the system-bath 
interaction is strong. We analyze the analogous behavior 
in the non-equilibrium spin-fcrmion model |(5J) in section 
V. 



where v and a n (n = L, R) are dimensionless system- 
bath coupling strengths 

v = np(e F )V L ,R; ai = np(e F )V L ,L; ct 2 = np(e F )V RiR , 

(6) 

and p(e F ) is the density of states at the Fermi energy 
of the L, R reservoirs. Out-of-equilibrium, Au 7^ 0, the 
phase shifts are complex numbers given by [36[ , 



tan 8 j 



tan#R = a 2 — % 



1 - ia 2 ' 



1 + iai 



(7) 



Since the reservoirs density of states weakly varies around 
the Fermi energy, Di > A/i, the phase shifts are approx- 
imately energy independent, and are all calculated at the 
Fermi energy cf [36j . For simplicity, throughout the pa- 
per we typically consider the case of a = a± = a 2 <C v, 
and take V n>n > to be real. We note however that our main 
results, in particular the appearance of Marcus-type be- 
havior in the non-equilibrium regime at strong coupling, 
can be rederived using other variants of this model sys- 
tem with no limitations on the strength of the diagonal 
V n ,n interactions, as well as for the spin-resonant-level 
model of Ref. see Appendix A. 

Under these simplifications, the non-equilibrium phase 
shifts are given by 



tan Sl = v 2 (i — a), 
tan<5ft = — v 2 (i + a). 



(8) 



For a — the inverse tangent in Eq. ([5]) has a branch 
cut, conventionally placed at (— ioo, — i] and [i, ioo). For 
this special case, 



(9) 



The weak potential limit therefore corresponds to 5 + = 
— 5- ~ v and Sl = —5 R ~ iv 2 , so \5l,r\ "C \5±\. How- 
ever, as v — * 1, 5l,r diverges, whereas the equilibrium 
phases 5± are finite. 



III. NONADIABATIC DYNAMICS 
A. Overview 

We are interested in the reduced density matrix pd(t) 
in the space of d occupancy. This is defined in terms of 
time evolution from an initial condition pi at time t = 0, 



Pd(t) = Tr leads [e 



- iHt p z e mt ] 



(10) 



If the parameter A in Eq. ^ vanishes, the problem 
is just electrons in a time-independent potential, and a 
closed form analytical solution exists. If A 7^ 0, pd{t) 
may be expressed as an expansion in A. Evaluation 
of any term in the expansion entails solving a problem 
of electrons in a time- dependent field. In equilibrium, 
an essentially exact closed-form solution exists, and the 
main problem is to re-sum the series in A [4|| |4t| . For 
non-equilibrium problems an analytical expression is not 
known. In this paper we present a detailed numerical 
evaluation of some low order terms in the expansion for 
A. The essential features are revealed by the Golden Rule 
decay rate, obtained by assuming that (i) = at 

time t — 0, and (ii) pi is the density matrix corresponding 
to the ground state of H with nd{+) = 0, and (iii) that 
an expansion to C(A 2 ) suffices. This level of description 
is equivalent to the " non-interaction blip approximation" 
(NIBA) in the standard spin-boson model [l|, [D, l5fj and 
yields the (nonadiabatic) Fermi Golden Rule for the for- 
ward (+) and backward (—) transition rates between the 
spin levels as [13, HH, H3, [53| 



-j 25R ^ e± iBt C f (t)dt; 



(11) 



where T denotes time ordering, H^{t) 
e %H B tjj(ft e -iH B anc j j s g roun d state energy 
of the two uncoupled reservoirs. The Hamiltonians H± 
are defined in Eq. ([3]), refers to the real part of the 
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integral, and the trace is performed over the electronic 
degrees of freedom. For convenience, the term including 
the energy bias B is taken outside of the trace. 

The object of our calculation is therefore the correla- 
tion function C/(t), which should be evaluated for non- 
equilibrium conditions covering time scales from Dt ~ 1 
up to Apt 3> 1. D, an energy of the order of the Fermi 
seas bandwidth, and the potential drop A/u, specify two 
inverse time scales in the problem, where we typically 
work in the limit of Ap -C D. Note that, unlike the 
equilibrium case, there is no exact analytical approach 
to calculate C/(i) valid for all time scales. Approximate 
analytical approaches and exact numerics may be per- 
formed, as discussed below. 



of the scattering matrix, was later generalized to include 
finite temperature effects 4l|. An exact formal determi- 
nant solution was presented in Ref. [42j for the study of 
tunneling in a non-equilibrium electron gas. 

A formal solution for C/(t) is obtained from the linked 
cluster theorem (valid also for non-equilibrium problems) 

MM, 



$f(t) = - [ dATr]/(i)£ A (M), 
Jo 



(13) 



with G x (t, t') the matrix Green function in the space of 
the leads for Eq. (TT|), but with Hsb — > A-ffss- For this 
model G solves the Dyson equation 



B. Short and long-time asymptotics; 
Non-equilibrium x-ray edge problem 

The correlation function Cf(t) [Eq. (fTTj) ] is a crucial 
element in the theory of the x-ray edge problem, an ef- 
fect originating from the many body response of a Fermi 
system to the fast switching of a scattering potential, e.g. 
the creation of a core hole [32, 33]. The x-ray edge Hamil- 
tonian is a simplified version of the spin-fermion model, 
Eq. @, with a static subsystem that is either empty or 
populated, 



H 



H 



(/) 



dU, 



Hsb — ^ Vk,n:k',n'al. n ak',n'd^d, 



Hf,, 



k,n=L,R 

Ho + Hsb- 



(12) 



Here dJ, d are creation and destruction operators of the 
core electron, and a\ n (<Xfc, n ) creates (destroys) an elec- 
tron in the n-th lead with momentum fc. The single band 
x-ray singularity problem was originally solved exactly in 
the asymptotic limit by Nozieres and De Dominicis (ND) 

In the last ten years there has been a growing inter- 
est in understanding the x-ray edge effect in the meso- 
scop ic regime I34l 13511 and for non- equilibrium systems 
[la M Ea Ella, SU, El, where the core hole couples 
to more than one Fermi sea at different chemical poten- 
tials. Standard equilibrium techniques, e.g. bosonization 
[251I261, 27], cannot be simply generalized to handle these 
non-equilibrium systems (see Appendix C). The first to 
address the non-equilibrium problem was Ng, who gener- 
alized the Nozieres-De Dominicis solution to include more 
than one Fermi sea with different chemical potentials 
[3^ ]. Ng demonstrated that the edge singularity could 
be described by generalized phase shifts which are real 
for equilibrium systems and complex when the system 
is driven out-of-equilibrium. Physically, complex phase 
shifts reflect the finite lifetime of a non-equilibrium sys- 
tem. More recently, Muzykantskii et al. [39l |40T ] for- 
mally solved the out-of-equilibrium problem using the 
Riemann-Hilbert approach. The result, given in terms 



gt(ti,t 2 )=g(ti,t 2 ) + J dT£(ti,T)Z A (r)g A (T,ta),(14) 

with V_ x (t) = V n , n >rid{+)(T) and the unperturbed 
Green's functions 



U - to 



(15) 



p is the reservoir density of states, taken to be the same 
for the L and R leads. 

In equilibrium, ND showed [32j that this equation can 
be solved exactly and the coupling constant integral per- 
formed, leading to 



c f (t) ~ (Dty 



(16) 



with D an energy of the order of the Fermi sea band- 
width, [3 = {5\ + S 2 _ )/7T 2 , and 5± defined in Eq. ©. For 
the a = model studied explicitly 



(3 



^ ( atan(z/) 

V 7T 



(17) 



Eq. (|16p also holds for the non-equilibrium problem at 
times Apt <C 1. 

At long times. Apt ^S> 1, the equation was solved by 
Ng [36| ; see also [HI . The coupling constant integral may 
similarly be performed, leading to 



C f (t) ~ e- TA ^(Apty, 



(18) 



withT = \S'l-S'^\/2ir, 7 



by Eq. 



Here 6'' 



-(^i + ^l)/^ 2 and S l,r given 
L, R) refers to the imaginary 



part of the phase shift. For the model studied numerically 
(a = 0) we have 



r 



1 

2^ 



In 



and 



7 = 



2tt 2 



In'' 



l-i/ 2 



(19) 



(20) 
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C. Intermediate time 

Although the long time and short time behav- 
ior is known essentially exactly, a transparent non- 
perturbative analytical expression for Cf(t) that encom- 
passes all time scales Apt and coupling strengths pV has 
not been developed. Indeed, in this work we argue that at 
strong coupling (y — ► 1) a different functional form dom- 
inates at intermediate times Apt ~ 1 — 10 where a promi- 
nent Gaussian decay emerges, Cf(t) ~ e -K(A^t) (£)£)-/? 
We first offer a perturbative calculation which suggests 
this result, and then present exact numerical simula- 
tions which prove this behavior. The dominance of the 
Gaussian behavior at intermediate times translates into a 
Marcus-type rate in frequency domain, with bias voltage 
activation (see discussion in section V), instead of tem- 
perature activation, as in the classical Marcus rate (see 
Appendix B). 

The correlation function Cf(t) can be evaluated us- 
ing the cumulant expansion 44]. Note that unlike the 
bosonic case, all cumulants contribute, 



CV(t) = exp 



n=l 



K n {t) = 
{-i) n rt 



o Jo 



dh / dt 2 ... / dt n (TF(t 1 )F(t 2 )...F{t n )) c , 



(21) 



where T denotes time ordering, F = 
a \n a k',n'-, and (...) c denotes a cumu- 
lant average. The first cumulant yields an energy shift, 
while the second term is given explicitly by (Dt > 1, 
a = Q) 



K 2 (t) = - 



\ j dh f dt 2 (TF(t 1 )F{t 2 ))c 
z Jo Jo 

2v 2 

—z- ln(Dt) 
n z 

1 — cos(A/ii) 
Apt 



2%Apt 



v 



Si(Apt) 



2 — [ 7e +ln(A M i)-Ci(A M <)]. 

7T Z 



(22) 



For details see Appendix D. The sine and cosine integrals 
are defined as Si(x) = f£ sl "^ dt, Ci(x) — j e + \n(x) + 

Jo C ° s( t K1 ^ and 7e = 0.5772 is the Euler-Mascheroni 
constant. 

This expression reproduces the weak coupling limits of 
the analytical results Eqs. (|T6|) and (fl~8|) at short and long 
times respectively, and provides an interpolation between 
the two times. In particular, the second line describes 
how the long-time dissipation TApt term is "turned on" 
as Apt increases from a small value to values much 
greater than unity. The first and last terms describe how 



the equilibrium orthogonality is turned off as Apt in- 
creases: [7 e + \n(Apt) — Ci(Apt) — \n(Dt)] is a function 
which interpolates between ln(i) for l/D < t < l/Ap, 
and a constant at Apt ^> 1. Note that in this model 
the leading logarithmic term at long times is ~ v A : con- 
sistent with the cancellation of logarithmic terms in the 
long time limit of Eq. (f2"2"l) , i.e. with the absence of a 
term which "turns on" the non-equilibrium power law. 
This cancellation does not necessarily occur at order v 2 
in other models, e.g. the spin-resonant-level model, sec 
Appendix A. 

We can clearly distinguish between three regimes in 
Eq. (EH): 



Cf(t) 



-(VA/ji) 2 /2tt 2 x £-2i/ 2 /tt 2 

-l/ 2 A/it/7T 



Apt < 1 

Apt ~ 1 (23) 

Apt > 1. 



While the first (equilibrium) limit and the third regime 
are well established in the literature [lj, [36|, HH, |4(| , the 
intermediate domain, leading to an interesting new dy- 
namic has not been discussed. In the strong coupling 
limit the Gaussian behavior may have a dominant ef- 
fect on the relaxation, as discussed below. We would like 
therefore to phenomenologically extend the second cumu- 
lant expression, Eq. (f2"2"|) . to larger phase shifts (strong 
coupling). 

Perturbative expressions analogous to Eq. (f2"2")l moti- 
vated Mitra and Millis [l(| to propose an interpolation 
function constructed by replacing the factors of v 2 in the 
expression above by the exact phase shifts. For the model 
considered here their procedure leads to 



*/(*) = 

, \&'L-&'k\ 



(6 2 + + S 2 _ 



(Apt) 



Si(A/rt) 



1 — cos(A/it) 



Apt 



[ln(l + iDt) - 7 e + Gi(Apt) - ln(A/xi)] 



<() 2 4- i) 2 ) 

L V^[7e-Ci(A M t)+ln(A M t)]. 



(24) 



Note that our approximation for the scattering poten- 
tials, ol\ — a 2 ~ 0, implies that there is no Fumi en- 
ergy shift. At the short time limit, Apt <C 1, the factor 
[\n(Apt) + 7e — Ci(A/^t)] dies out, leading to the correct 
equilibrium behavior (|16[) . In contrast, at long times the 
cosine integral diminishes, which implies that the dynam- 
ics is ruled by an exponential decay with a rate constant 
ApT, [Eq. (|19p ], modified by a power law term i 7 , Eq. 

Our numerical results, to be presented below, show 
that at weak to moderate coupling, v < 0.5, the correla- 
tion function and the resulting transition rates are well 
described by expression (|24|) . However, Eq. (|24[) is found 
to be a poor approximation at strong coupling. Instead, 
at intermediate times Apt ~ 1 we return to Eq. (|22p and 
replace the weak coupling phase shift by the equilibrium 
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strong coupling phase shift, v — > atan(z/). The physi- 
cal picture is that on these time scales the phase shifts 
are essentially still the equilibrium ones. Only at longer 
times A/it 1 the non-equilibrium dynamics is reflected 
in the complex phase shifts ([8]). This conjecture yields 

^(A/zt ~ 1) = $ eg (t) + $ neg (t) + iE s t, (25) 

where the equilibrium function is the same as in the zero 
bias case [321 ] . 

$e,(t) = /31n(l + iDt); 

while the non-equilibrium term provides a quadratic time 
decay 

= «(A/it) 2 ; 

k = (^+5) = ^M. (27) 



4tt 2 



2tt z 



Notice that the prefactor k depends only on the scat- 
tering potential V n , n '- The last element in Eq. ([25]) is 
the energy shift E s . We assume that it is given by the 
equilibrium limit of the Fumi's theorem, 



E s = -[6. 



(28) 



For a = the energy shift is zero. 

Similarly to expressions [j2"5 ]) -([25 ]) . Eq. (gj) gives the 
first correction to the equilibrium result which is propor- 
tional to (A/it) 2 , but in contrast to these equations, the 
coefficient involves the non- equilibrium exponent, and in 
fact does not provide a wide regime of t 2 behavior. Our 
numerical simulations, presented below, support expres- 
sions ([25" ]l -([2"? ]) in the broad window A/it ~ 1 - 10. We 
have not been successful in constructing a general ana- 
lytical expression, valid on all time scales and coupling 
strengths. It is possible that consideration of the fourth 
cumulant may yield some insight here. 



IV. NUMERICS 

A. Methods 

The fermionic correlation function C/(t) can be di- 
rectly calculated by expressing the zero temperature 
many body average as a determinant of the single particle 
correlation functions 144, y54j 



G f (t) 



= det[(pk,n;k',n>{t)] k<k n. k , <k n> J 

= (fc^le^+V^lfcV). (29) 

Here H± = ^ h±, where h± are the single particle 
Hamiltonians for the individual conduction electrons. 



\k, n) are the single particle eigenstates of Hg, and the 
determinant is evaluated over the occupied states, is 
the Fermi energy of the n-th reservoir. In our numer- 
ical calculations we have used a Lorentzian density of 
states, with tails that are long enough to eliminate ar- 
tificial reflections from the boundaries. The Lorentzian 
function is centered around the equilibrium Fermi energy 
with a full width at half maximum Dl, Eq. Q. This 
quantity sets energy and time scales in our simulations. 
We have typically used Dl = 4.5 for the two reservoirs, 
A/j,/Dl < 0.1 and v = 0.1 — 1. We also take the diago- 
nal coupling to be zero (a = 0) in all of our simulations, 
unless otherwise stated. For these parameters, we have 
found that for short-time evolution (A/it < 15), even for 
strong coupling, it is satisfactory to model the fermionic 
reservoirs using ^400 states per bath, where bias is ap- 
plied by depopulating one of the reservoirs with respect 
to the other. 

We can also employ the renormalization group (RG) 
method, originally developed by Wilson for the calcu- 
lation of the thermodynamic properties of the Kondo 
problem [551 ] . for the numerical solution of the non- 
equilibrium x-ray edge problem. In equilibrium, Oliveira 
et al. [56|, H3] have used the RG technique to calculate 
the x-ray absorption spectrum. This study can be gen- 
eralized to include two reservoirs with different chemical 
potentials by following a three-step procedure: (i) de- 
fine the conduction bands on a logarithmic scale, (ii) 
convert the (isolated) reservoir Hamiltonians into semi- 
infinite tight binding chains, as is done in Refs. [56l l57l] . 
In this representation the impurity couples the chains' 
first levels, (iii) build the Hamiltonians H± in the new 
basis, first including the occupied levels of the L and R 
reservoirs, then adding the empty levels. The determi- 
nant (|29[) is performed over occupied levels only. 

In equilibrium, the RG technique is highly advan- 
tageous over constant/Lorentzian discretization meth- 
ods, as it converges rapidly to the continuum limit even 
for gross discretization. For small voltage differences 
(Afi/D < 10~ 2 ) this method nicely reveals the crossover 
of Cf(t) from equilibrium to non-equilibrium behavior 
with increasing bias. In contrast, for large bias the 
Lorentzian discretization is more convenient, since en- 
ergies far from the Fermi energy are not well represented 
within the RG technique. Wc present a numerical exam- 
ple in Fig. [TJ demonstrating the strength of the RG ap- 
proach over standard linear discretization for systems in 
equilibrium. The RG technique provides stable dynam- 
ics for long times (full line), where constant discretization 
fails (dotted line), yielding an artificial rise of the correla- 
tion function due to discretization errors. The theoretical 
value of (3 = 2atan 2 (^)/7r 2 = 0.0020 for v = 0.1, nicely 
agrees with the numerical slope of 0.0019. Deviations 
are due to the sharp energy cutoff used at D —l, with 
the conduction band energies extending from —Dq to Dq. 
We also present the results of an RG calculation with a 
very small voltage drop (dashed line), where linear dis- 
cretization would require a very fine grid. 
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In this work we typically focus on systems far from 
equilibrium, A^i/D ~ 0.1. Since the RG method sam- 
ples the Fermi sea states predominantly near the Fermi 
energy, while high energy states are under-represented, 
we find the Lorentzian discretization to be more conve- 
nient. 



B. Results: C f (t) 

Representative results are displayed in figure [2] The 
main plot presents the logarithm of the correlation func- 
tion |C/(t)| at strong coupling v = 0.95 for an applied 
voltage Afj, — 0.24. Three different regimes are clearly 
identified: a power law decay at short times A/it < 1, 
see lower left inset (a), an exponential decay at long 
times A /it 3> 1, and remarkably, an intermediate regime 
1 < A/xt < 10 of approximately Gaussian behavior [up- 
per right inset (b)]. The short and long time behaviors 
are consistent with the theoretical results. The interme- 
diate time quasi-Gaussian regime is a new finding with 
important consequences. We analyze the short time dy- 
namics A/it < 1, enlarged in Fig. HJa), by fitting the data 
to the analytic expression In C/(t) ~ — j3\sx{Dt), This 
provides an effective bandwidth D = 6 and a decay con- 
stant [3 = 0.13 consistent, within numerical errors, with 
the theoretically expected (3 = 2[atan(;/)/7r] 2 « 0.12. We 
can also fit the intermediate time behavior, shown in Fig. 
Hfb) , by a Gaussian function In C/(t) ~ — «;A/i 2 t 2 which 
yields the prefactor k=0.03. 

Figure [3] presents a more detailed examination of the 
Gaussian behavior, showing that at both short and in- 
termediate times, the data can be well described by the 
approximate function 

y G (t) = (Dt)-^e-< A ^ 2 , (30) 

with j3 the theoretically predicted short time (equilib- 
rium) exponent (6+ + <?_ ) /ir 2 . The inset proves that the 
data follows the same linear trend when plotted as a func- 
tion of (A/it) 2 , with a slope of k ~ 0.03. This value 
nicely agrees with the constant predicted by Eq. (f2"T)) . 
k = atan(;/) 2 /27r 2 =0.029 [y = 0.95). 

Fig. [4] provides more insight by deconstructing the ob- 
served time decay of C/(t) into the equilibrium power 
law and non-equilibrium Gaussian components. Another 
important observation deduced from Figs. [2 [3] and [4] 
is that the correlation function decays to ~ 0.1 its ini- 
tial value by the time the exponential decay begins to 
dominate. This implies that the Gaussian behavior gov- 
erns the rate constant at strong enough coupling, lead- 
ing to a voltage activated regime analogous to the high- 
temperature semiclassical polaron transport regime. We 
call this "fermionic Marcus" behavior. 

Fig. [5] presents the evolution of the correlation func- 
tion C/(t) as coupling strength v is varied from weak to 
strong. All other parameters are the same as in Fig. [5] 
For all coupling strengths the short time logarithmic and 








FIG. 1: Decay of C/(t) as computed by both the RG ap- 
proach and standard linear discretization for v — 0.1. Con- 
stant discretization with 200 states per band, A/i = (dot- 
ted), showing an artificial rise of the correlation function 
due to discretization errors; logarithmic discretization (RG) 
with 100 states per band, A = 1.1, A/i = (full); logarith- 
mic discretization (RG) with 100 states per band, A = 1.1, 
A/i = 0.007 (dashed). A is a logarithmic scale parameter, 
where for each conduction energy there are states with ener- 
gies e/A m , m = 1,2.... Inset: The finite bias case exhibits an 
exponential decay at long times. The slope agrees with the 
theoretical value, Eq. (I19p . The Fermi sea bandwidth is 2Dq 
in all plots with Do = 1. 



the approximate long time exponential behavior are ob- 
served. However, as the coupling strength is increased, 
increasingly wide intermediate regime is observed. We 
have verified, by an analysis similar to that shown in the 
lower left inset of Fig. [51 that the short time behavior is 
always a power law with the theoretically predicted ex- 
ponent (3 = 2[atan(z/)/7r] 2 . Also note that while at weak 
coupling (y < 0.5) the correlation function weakly de- 
cays before the turnover to an exponential decay takes 
place, for very strong coupling, v > 0.9, the dynamics is 
critically controlled by the Gaussian form, as the correla- 
tion function has decayed to zero before the exponential 
decay takes place. This implies that the resulting decay 
rate [Eq. I|ll[) ] essentially shows different characteristics 
in these two regimes. 

We now systematically explore the Gaussian decay at 
intermediate times, and the exponential decay rate at 
long times, and compare the numerical coefficients T and 
K with the theoretical values, Eqs. (fTTj)) and (|2"T[) , re- 
spectively. This is done by calculating the correlation 
function C/(t) for coupling strengths v = 0.1 — 0.95 (see 
Fig. [5J, then extracting both the quadratic intermediate 
slope kA/j, 2 and the long time exponential slope TA/i. 
Fig. [5] presents these coefficients showing excellent agree- 
ment with the values predicted from the phenomenolog- 
ical ansatz, Eqs. (|2l) ]) -(|27 ) . 

Next, in Fig. [7] we examine the crossover to the an- 
alytic long time behavior, Eq. (|24p . We compare the 
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FIG. 2: The correlation function C/(t) for Art=0.24, ^=0.95, 
manifesting a power law decay at short times Afit < 1 (a) 
(notice the log- log scale), a Gaussian decay at intermediate 
times A/j,t ~ 1 — 10 (b), and an exponential decay at long 
times A/j,t 3> 1. 
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FIG. 4: Time dependence of the correlation function |C/(t)| 
for Att = 0.24, v = 0.95. The fitting function y G = 
e - K (A M t) 2 p^-/3 ( dotted ) with K = q Q3 ) D = 6 and /3 = 0.13 
compared to exact numerical solution (full). The Gaussian 
(dashed-dotted) and the power law (dashed) parts of j/g are 
also displayed separately. 




FIG. 3: Evidence for the Gaussian decay at intermediate 
times A/j,t ~ 1-10 and strong coupling v = 0.95, A/j, = 0.16 
(full), All = 0.24 (dotted) and Afi = 0.32 (dashed-dotted). 
The inset, which includes the three lines one on top of 



the other, reveals that Cf{t)v 
k ~ 0.03. 



-K,(Anty 



, (3=0.13, with 



numerical correlation function with two functions: the 
approximate fitting function yc defined abovej_and the 
long time perturbation theory result of Ref. given 
by exponentiating Eq. We refer to this second 

function as tje- We see that he describes the data well 
at long times, but that as the coupling strength is in- 
creased, the range over which the Gaussian description 
applies increases. This feature can be qualitatively de- 



scribed by the approximate crossover function 



c { ; pp \t) 



exp ■ 



4 1 1/2 



4k 2 



2k 



(31) 



which captures the crossover from a Gaussian dynam- 
ics to an exponential decay. An increase of v leads to 
a strong enhancement of T, while k reaches saturation, 
resulting in a counterintuitive lengthening of the range 
of the intermediate Gaussian dynamics with increased T. 

In summary, we have shown that the crossover between 
equilibrium (A/xf <C 1) and non-equilibrium (A/it ^> 1) 
behavior is described by a regime of Gaussian relaxation 
negligible for weak coupling, but for strong coupling ex- 
tending over the wide range 1 < Afi < 10, with parame- 
ters determined by the equilibrium exponents. In Section 



9 



0.4- 
0.2- 



0.5 1 1.5 2 2.5 

ln[(1+v 2 )/(1-v 2 )] 



3 3.5 



0.03 
0.02 
0.01 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 



FIG. 6: Testing the validity of Eqs. (JT9J) and J27)) for de- 
scribing the long time and intermediate time behavior, re- 
spectively. Top: The relaxation rate V. Numerical results, 
calculated from the slope of \n[Cj{t)t^] vs. A[it at long times 
(squares); analytical results using Eq. (|19[) (dashed line). 
Bottom: The coefficient k. The numerical slope of ln[C / (i)t ] 
vs. A(i 2 t 2 at intermediate times (squares); analytical results 
using Eq. (|27[l (dashed line). These data were computed with 
A/x=0.24, and (5 was extracted from the short time dynamics 
for each value of v. 



V we examine the consequences for spin relaxation. 
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C. Orthogonality and anti-orthogonality 

We focus next on the power law contribution to Eq. 
(fT8| . Unlike the standard equilibrium case, where the 
system always experiences dephasing, 7 < |58(, in our 
model the power law term in Eq. (|18[) acquires a posi- 
tive exponent 7 > 0, enhancing the correlation function, 
see Eq. (j2"0")l . We refer to this situation as an "anti- 
orthogonality" effect. For a general system-bath coupling 
model, Eq. (f24|) reveals that 

C / (A^»l)^e- rA ^; i = -{5l+5l)/n\ (32) 

with complex, non-equilibrium, phase shifts given by Eq. 
d?]) (3(| ■ It is clear that in the special limit of zero diago- 
nal interactions (a = 0), the phase shifts 8 l ,r are purely 
imaginary and (5 L + 6%) < for all values of v, leading 
to 7 > . In contrast, for large diagonal coupling we 
typically find that 7 < 0, which is the standard orthogo- 
nality behavior. The anti- orthogonality effect is therefore 
a footprint of a non- equilibrium situation. 

We next turn to a numerically exact exploration of the 
anti-orthogonality effect. Fig. [8] shows the correlation 
function Cf(t) at long times A/i£ 3> 1 when expression 
(|24| holds. We numerically extract the long time slope 
TA/x, and recover the weak power law dependence by 
multiplying the correlation function by the inverse of the 
exponential decay. The standard orthogonality effect is 



FIG. 7: The turnover from a Gaussian decay to an expo- 
nential relaxation. Top: Comparison between the numeri- 
cal correlation function Cj(t) (full) and the fitting functions 
vg (dotted) and ve (dashed) defined in the text, v — 0.1, 
A/x=0.24. Inset: Exposing the turnover by taking away the 
short time and intermediate time terms t _l3 e~ ^ 4 with 
,5=0.002, k = 5 x 10~ 4 . Bottom: Same with v = 0.95, leading 
to /3 = 0.13, k — 0.03. While ya explicitly includes the Gaus- 
sian decay, correct to A/it ~ 10, the function yE captures the 
correct slope at longer times. 



presented in panels (a)-(b) for a — 0.5 and v = 0.3, for 
which 7 = —0.038. When a = and v — 0.8 the anti- 
orthogonality effect clearly manifests itself with 7 = 0.12 
(c)-(d). Interestingly, the correlation function at long 
times shows a complicated behavior, more complex than 
that predicted in Eq. as evidenced by the mild 

deviations from strict power law behavior displayed in 
Fig. ETd). 

Fig. [5] presents an " orthogonality-anti-orthogonality" 
map as a function of the diagonal (a) and nondiago- 
nal (v) couplings using the general expressions of Eq. 
([7]) with a = oti = a 2 - We find that for large a, 
— 1/2 < 7 < 0, manifesting the standard orthogonal- 
ity effect. For large nondiagonal interactions typically 
anti-orthogonality may be observed. 

In addition to the long time exponential decay and 
anti-orthogonality behavior, non-equilibrium dynamics 
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FIG. 8: Orthogonality and anti-orthogonality effects in the 
non-equilibrium system A/i = 0.4. (a)-(b) ^=0.3, a=0.5, 
manifesting the standard orthogonality effect (7 < 0). (c)-(d) 
v=0.S, a—Q, revealing the anti-orthogonality effect (7 > 0). 
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FIG. 9: Map of orthogonality [(S 2 L + 5 2 R ) > 0] and anti- 
orthogonality [(Si, + 8 R ) < 0] behavior using expression (J7J) 
with a = ot\ — ct2- 
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FIG. 10: Resolving the complex part of the power law ex- 
ponent 7. Qi = 0, a2=0 (full); ai = 0.1, «2=0 (dashed); 
on = 0.2, a 2 =0 (dashed-dotted); an = 0.2, a 2 =0.2 (dotted); 
Afi = 0.24 and v = 0.95 in all plots. 



reveal an additional decaying contribution which is ex- 
pected to oscillate at longer times. We did not suc- 
ceed in fitting I(t) to the stretched-oscillatory function 
e *7 lnt^ indicating that at strong coupling the dynamics 
is more involved. Finally, we note that, consistent with 
our observations above, equilibrium effects dominate up 
to A/it ~ 10. Only at longer times I(t) begins to deviate 
from unity due to the emerging influence of the imaginary 
term 7". 

Though the imaginary term 7" can strongly affect the 
correlation function Cf(t), (Fig. [T0|) . its practical con- 
tribution to the Golden Rule rate is small. We find that 
for weak to intermediate coupling, 7 <C 1, leading to 
I(t) ~ 1. On the other hand, for strong coupling, 7" 
manifests itself only at long times A [it > 10, when the 
correlation function has essentially decayed to zero. 



may be reflected in the appearance of complex power law 
exponents [36j . When the spin impurity is symmetrically 
coupled to the two leads (ai =02), the phase shifts are 
complex conjugates, see Eq. (J7J, and 7 is always real. 
This situation was discussed above. In contrast, asym- 
metric systems may acquire a complex coefficient with 
7 = 7' + vy", a direct outcome of a non-equilibrium sit- 
uation. 

The imaginary contribution to 7 is resolved in Fig. 
I10I Since at weak coupling the imaginary term 7" is 
very small, we investigate a strong coupling system with 
v = 0.95. Motivated by Eq. (|24|) . we assume the 
generic form C f (t) = \C } (t)\e ltt e^" ln *. We numeri- 
cally extract the phase factor e, then plot the func- 
tion I(t) = Cf(t)/\Cf(t)\e~ %et for different diagonal cou- 
pling strengths. As expected, in symmetric situations, 
I(t) ~ 1. In contrast, asymmetric systems (a% 7^ 0,%) 



V. RELAXATION 

A. Qualitative discussion 

In this section we present calculations of the nonadi- 
abatic relaxation rates Tf(oS) defined in Eq. (TTTj) . The 
given physical model corresponds to two states separated 
by an energy uj which depends on the bare level split- 
ting B and on renormalizations arising from the cou- 
pling to the leads, (u < 0) corresponds to the up- 
scattering rate describing transitions from the lower level 
to the upper level, while r^(w > 0) corresponds to down- 
scattering. In equilibrium at T = 0, (lo < 0)=0, i.e. 
there is no up-scattering. At temperature T > 0, the 
detailed balance relation of equilibrium thermodynam- 
ics implies (~uj) /T^ (uj) = e~ w / T . In this section we 
examine the rates in the non-equilibrium situation. We 
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show that the Gaussian form of the correlation function 
Cf(t) which occurs at strong coupling has important con- 
sequences for the physics. 

Before discussing our results in detail, we establish the 
relevant energy scales. The general expression, Eq. (jlip . 
may be written (neglecting overall factors) as 



r 

r+M = K/ 

Jo 



dt 



(iDt) 



iwt—A(i) 



(33) 



Here f3 e f / (t) is an effective exponent which changes from 
the equilibrium power (3, Eq. (JTTJ) , to the non-equilibrium 
value (3 neq [defined as (—7) in Eq. (|2"0)) ]. as Afit changes 
from less than unity to much great then unity, uj is the 
physical energy level difference, given by the sum of B 
[Eq. @] and the level shift arising from the system-bath 
coupling, and D is an energy scale of the order of the 
Fermi sea bandwidth. 

The naive assumption [Icl ] is that the only important 
energy scale is the relaxation rate given by the current 
flow across the system, TA/i. In fact the numerical and 
analytical results presented in the previous sections in- 
dicate that the situation is more subtle. At short times, 
A(t) ps re(A/ii) 2 whereas at long times A(t) — > TA/ii. 
The interplay of re, which is proportional to coupling 
strength v 2 at weak coupling but saturates at strong cou- 
pling [Eq. (|27|) ]. and T which is proportional to v 2 at 
weak coupling but diverges at strong coupling [Eq. Q19p], 
gives a richer behavior. 

Appendix E gives details of an asymptotic analysis of 
Eq. (f3"3"| . This analysis reveals that to discuss the relax- 
ation rate one should distinguish strong and weak cou- 
pling. In the weak coupling limit, there are two relevant 
scales, YA\x and A/i (the latter multiplied by various 
factors which are in practice fairly close to unity). For 
uj > Afj, we get the equilibrium down-scattering rate; for 
\uj\ < Afj we find a nontrivial approximately Lorentzian 
behavior, and for uj < —Afj we reproduce the e _s ^ ln ^ 
behavior found by Mitra et al. [lOj]. Specifically, 



= r 



TAfi Afj 
mT 



T 2 Afj 2 

■ ln(w/A/i). 



P 

; ui > A/i 

< uj < A/i 
Ait 



«<-j^- (34) 



be important: i/kAli and TA/i. We find 

= ni-Pf^(%Y-, u>V&» (35) 



UJ 

f[(l-/3)/2]cos^ (VnAjiV, 



e 4k(Ah) 



2y/K,AfJ, \ D J 

H< V 
T cos(7r/3) ( 2reA/i 



uj\ < y/nAfi (36) 

2\P 



2y/nAn \ Dlu 

-TAfi < uj < -V^A^ (37) 



e a ^ 



ln(td/A/i). 



uj < -TAfi 



(38) 



The Gaussian behavior found at intermediate frequency 
scales is a consequence of the wide regime of t 2 behav- 
ior found in the time evolution function, and may be 
roughly understood as the Fourier transform of e -K ( A/it ) 
although as the results of Appendix E show, this argu- 
ment must be treated with some care. 

We call Eqs. p6 ]) - ([3"7) the "fermionic Marcus rate," 
the analogue of the classical Marcus result for spin-boson 
systems (|B12I) , which holds in non-equilibrium situations 
at strong coupling. This expression indicates that the 
voltage activates the absorption rate, similarly to the role 
of temperature in the bosonic case. The result differs 
from the bosonic solution (Appendix B) in some impor- 
tant aspects: (i) In the fermionic case the Gaussian decay 
is modified by a weak power law term, (ii) For bosonic 
systems the activation factor depends on the tempera- 
ture as InTf oc T^ 1 , while for fermionic systems we get 
a voltage squared activation, lnT/ oc A[i~ 2 . Therefore, 
there is no simple linear mapping between temperature 
and voltage drop in the strong coupling regime. We note 
however that the classical Marcus rate is applicable in 
the high temperature limit (see Appendix B), while we 
typically assume here that Afj, <gC D. Therefore, in our 
system the energy window for reorganization processes is 
the bias voltage, rather than the full bandwidth D. Thus, 
we may interpret the kA/i 2 factor in the denominator of 
the Gaussian decay (|37[) as a reorganization energy of 
the non-equilibrium fermionic system, A/ ~ nAfi, mul- 
tiplied by the driving force Ff = A/i. In contrast, in 
the equilibrium spin-boson model, reorganization ener- 
gies are of order of the cutoff frequency, oc uj c , and 
the driving force for absorption processes is temperature 
Fi, = T. Qualitatively, both models then recast to the 
familiar Marcus- like form, T oc l XF 59]. Further, 
both the fermionic and the standard Marcus behaviors 
share similar qualitative features such as the existence of 
an inverted regime, as discussed in the next section. 



B. Numerics: rates, population and current 



Here f (x) is the complete Gamma function. Note that 
the formulae match at uj ~ A/i because in weak coupling 
(3 ~ T < 1 leading to T+(w - Afj,) oc 1/uj. In the 
strong coupling limit, two frequency scales turn out to 



We numerically evaluate the integral ([33]) using the 
coefficients r, 7, (3 and re as determined by the coupling 
strengths, Eqs. fTS]). PD|1 PS]) and ^ respectively. For 
convenience, we disregard the multiplicative factor A 2 /2. 
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FIG. 11: Golden Rule relaxation rate calculated for the weak 
coupling limit v = 0.5 (/3 = 0.043, F = 0.08), A/i = (dashed 
line); A/i = 0.24 (solid line); A/i = 0.48 (dotted line). The 
inset presents an expanded view of low frequency regime. 
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FIG. 13: Spin polarization as a function of potential bias A/i 
for different energy biases u = 0.1 (full), id = 0.2 (dashed), 
id = 0.3 (dashed-dotted) for n = 0.03, = 0.13 {v = 0.95). 
Inset: Mean field calculation of the current / <x (nti(+)) 2 A/i 
for the same frequencies as in the main plot. 
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FIG. 12: Fermi Golden Rule rate, Eq. (I33p . as a function of 
frequency, ft = 0.03 and j3 = 0.13 (^ = 0.95). Inset: Voltage 
activated excitation rate (u < 0). 



The main panel of Fig. I 111 shows on a semi-logarithmic 
scale the relaxation rate computed numerically for the 
relatively weak coupling v = 0.5 (/3 = 0.043, T = 0.08, 
non-equilibrium exponent (3 neq = —7 = 0.013) and two 
choices of chemical potential, A/i = 0.24 and A/i = 0.48. 
Also shown as the dashed line is the T — equilibrium 
result. The inset shows an expanded view of the small 
frequency regime, demonstrating the Lorentzian behav- 
ior. We clearly observe the three regimes as discussed in 
Eq. (f3"4")l : For small frequencies (large bias voltage) the 
spin levels are approximately degenerate, and the rates 
are symmetric around u> — (inset). In the opposite 



\ll>\ > A/i limit, the absorption rate is practically zero, 
while the emission rate approaches the equilibrium limit. 
In between, a voltage activated excitation behavior is re- 
vealed. 

We analyze next the strong coupling limit. Fig. [T2l 
shows that the excitation process is activated by a finite 
potential difference as prescribed by Eq. (|3T|) . More 
quantitatively, the inset verifies that the relationship 
lnFj (w < 0) oc — (ui 2 /Afj, 2 ) holds. Similar to classical 
bosonic Marcus rate [45j, an inverted regime appears for 
the fermionic system. However, in the present case the 
rate in the inverted regime decays weakly as a power law 
rather than as a Gaussian. At large frequencies, u 3> A/i, 
equilibrium behavior is observed where rjr (w < 0) ap- 
proached zero, and rj^(u > 0) becomes insensitive to 
voltage. We have also calculated the Golden Rule rate 
using the numerical correlation function (depicted e.g. 
in Figs. [Hand [2]), instead of the approximate analytical 
function in (f3"3"|) , and find that the results agree perfectly. 

We now turn to a study of the spin polarization. In the 
incoherent tunneling regime, for small tunneling parame- 
ter A, the populations of the two levels obey a Markovian 
balance equation 

P + =TJP_-T+P + ; P_+P+ = l, (39) 

with the absorption and emission rates given by Eq. (|33|) . 

The polarization (a z ) = P + — P_ = f + L , shown in Fig. 

1131 manifests a transition from a fully polarized system 
(a z ) ~ —1 to an unpolarized system (a z ) ~ as A/i 
is increased. Typically, we find that the crossover takes 
place when the energy bias u> becomes comparable to 
the bias voltage. While at high frequencies, \ui\ ^> A/i, 
rt(w < 0) ~ 0, leading to full polarization, at very large 
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pose the propagators, including all spin-flip events along 
the time ordered contour, and obtain, e.g. for the spin 
up population 0, El[ , 



c 4 (ti.t 2 ,t 3 ,t 4 ) [— i r~i , A4 
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FIG. 14: Schematic representation of spin-flip evens on the 
Keldysh contour. Plotted are examples for particular two, 
four and six spin-flip processes, respectively. 
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bias the emission and absorption rates are comparable, 
resulting in equal population of the two levels and zero 
polarization. The Gaussian activation term in Eq. (|3"T|) 
is therefore reflected in the enhancement of polarization 
with bias voltage. 

It is also interesting to note that the electron current 
through the system, calculated at the level of mean field 
theory, / oc A/i(nd(+)) 2 , is strongly suppressed for weak 
bias, Afj, < u>, see inset of Fig. (T3J In contrast, for very 
large bias, (rid(+)) = 1/2, and the current increases lin- 
early with A/i. Therefore, it is the intermediate regime 
of lo ~ Afi that manifests prominent nonlinear current- 
voltage characteristics, emerging due to the interplay be- 
tween the Gaussian relaxation and the power-law dynam- 
ics. We can compare our results to the weak coupling 
Bloch-type rate equations of Gurvitz et al. |i60j which 
yield (a z ) = at long enough times, independent of 
voltage drop and energy bias. In contrast, Fig. [T3] re- 
veals a rich dynamics in the strong coupling regime with 
a prominent dependence on system energetics and the 
non-equilibrium conditions. 



VI. 



BEYOND C(A 2 ): COULOMB GAS 
BEHAVIOR 



In this section we discuss a crucial ingredient of the 
physics of our model that allows for a description be- 
yond the Golden Rule [C(A 2 )] level. A formally exact 
solution for the impurity spin problem ([2]) can be writ- 
ten by a power series in the tunneling matrix clement A 
(6l| . Here, we restrict ourselves to an exact numerical 
investigation of the electronic correlation functions that 
appear in this power series to see if the usual " Coulomb 
gas" behavior is observed even when the system is out- 
of-equilibrium. In particular, the reduced density matrix 
of the spin impurity Pd(t) is given by 



Pd (t) = Tr e-^V(0)e 



AHtl 



(40) 



with forward and backward time evolution branches. 
Here p is the total density matrix, and the trace is per- 
formed over the reservoir electronic states. We decom- 



00 00 



iM*)i+) = EE(- 1 ) (fc+i) (?) 



» x 2k+2j 



x / dti 
Jo 



k=0 j=Q 

dt 2 k / dsi 
Jo Jo 



S 2 j-1 



ds 2 j 



xTr 



e iH + s 2 j e -iH-s 2 j 



x[e iH+t le -iH_t, 



^ill-si ^-iH+si 



(41) 



Here |±) are the up and down spin states, and H± is de- 
fined in Eq. ((4]). This expression was derived assuming 
that at t = the spin is in the pure state |+), and the 
(isolated) reservoirs are in their respective ground states 
(T = 0). It can be easily generalized to describe other 
initial conditions. Each time variable in Eq. (|4ip marks a 
particular spin-flip event. While the second order correla- 
tion function couples nearest neighbor events only, higher 
order correlations couple distant spin-flips, yielding a 
multiparticle interaction term. In equilibrium, this inter- 
action can be exactly written in terms of pair-wise con- 



tributions, Tr[...] tx exp 



_E fc <,(-l) fc+ ^/(ifc-^)j, Big- 
nificantly simplifying the computational problem. This 
is the celebrated Anderson- Yuval-Hamann (AYH) result, 
which leads to the interpretation of the Kondo problem as 
a one-dimensional Coulomb gas system [62] ]. In contrast, 
in the general non-equilibrium case, the exact structure 
of the interaction is not known for all times, and it is not 
clear whether higher order correlations can be exactly 
decomposed into pair- wise contributions fl3j| . 

Using the numerical technique discusses in section IV, 
we can exactly calculate, term by term, the correlation 
functions in Eq. (|4"Tj) . Specifically, we study three ex- 
amples of the processes of the order of A 2 , A 4 and A 6 , 
depicted schematically in Fig. Qj] 



C 2 {tiM) = ( e «*-*i e <ff+(*»-*Oe"M*-*0e- <ff -*) 
Ci(ti, t 2 , ts, ti) = 

/gUT-tl e iH+(t 2 -ti) e iH-(t 3 -t 2 ) e iH+(tt-t s ) e iH_(t-t 4 ) e ~iH-t\ 



Ce(ti,t 2 , £3, £4, t&, te) 



(JH-ti e iH+(t 2 -ti) e ii?_(t 3 -i 2 ) e iH+(U-ta) e iH_(i 5 -U) 
xe iff + (t 6 -t 5 ) e iH-(t-t B ) e -iH-t\ 



(42) 
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and compare the results to the Coulomb gas expressions 
C 2 {t u t 2 ) = C 2 (t 1 ,t 2 ) = C 2 (t 2 -t 1 ), 

C4,{tx...ti) — 

c 2 (t 2 - tjXhtu - h)c 2 (t3 - t 2 )c 2 (u - h) 
c 2 (h - h)c 2 (t 4 - 1 2 ) 

Ce(h...t 6 ) = C 4 (h..U) x 

C 2 (t 6 - h)C 2 {t 5 - t 2 )C 2 (t 6 - t 3 )C 2 (t 5 - U)C 2 {U - t 5 ) 
C 2 (t 5 - tx)C 2 (t 6 - t 2 )C 2 (t 5 - t 3 )C 2 (t e - u) 

(43) 

In particular, our calculations were performed assuming 
a regular interval r between spin flips. 

We have robustly checked that the AYH decomposition 
[62[ holds precisely at all times (greater than Dt > 1) 
and coupling strengths in equilibrium, as well as for all 
times out-of-equilibrium, for weak to intermediate cou- 
pling strengths, see Figs. [TMTGl Interestingly, the AYH 
decomposition breaks down for the strong coupling out- 
of-equilibrium situation, precisely in the time window 
where Cf(t) shows a broad Gaussian decay with time, 
as depicted in Fig. [5j Even in this regime, the pair-wise 
AYH decomposition holds asymptotically for long and 
short times. 

An important outcome of this observation is that the 
AYH Coulomb gas expression [63 |. which is exact in 
equilibrium, cannot be justified for intermediate times 
A/i ~ 1 — 10 for strong coupling to the leads in the out-of- 
equilibrium situation. This is because at strong coupling 
the effective short time behavior (which cannot be de- 
scribed by the Coulomb gas picture) practically extends 
to longer times of order A fit ~ 1 — 10. The Coulomb 
gas expression still holds for weak to intermediate cou- 
pling strength and at long times. This investigation lays 
the groundwork for an exact evaluation of the spin dy- 
namics via path integral techniques valid even when the 
Coulomb gas decomposition does not hold. This work 
will be reported in a future publication [461 ]. 
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FIG. 15: Testing the Coulomb gas picture for out-of- 
equilibrium situations. Comparison between the exact fourth 
order correlation function C4(ti,t2,t3,t4), Eq. (|42p (full), 
and the pair- wise approximation C4(ii, £2, i3j £4), Eq. (I43p 
(dashed), r = ti+i—ti is the distance between spin-flip events. 
All other parameters are the same as in Fig. [S] 
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FIG. 16: Testing the Coulomb gas picture for out-of- 
equilibrium situations. Comparison between the exact expres- 
sion sixth order correlation function Ce(ti,t2,t3,t4,,t^,te), 
Eq. (|42|l (full), and the pair- wise approximation 
C(j(ti,t2,t s ,t4,t 5 ,t(j), Eq. (gSJ (dashed), r = U+i — U is 
the distance between spin-flip events. All other parameters 
are the same as in Fig. [S] 



VII. SUMMARY 

In this paper we have undertaken a detailed study of 
the non-equilibrium dynamics of a small quantum system 
coupled to two electronic leads. This problem is of great 
interest for understanding dissipative effects in prospec- 
tive single molecules devices. The model studied here 
is generic enough to capture range of relevant relaxation 
motifs, while being simple enough for detailed investiga- 
tion. Our analysis combines analytical results with exact 
numerics, rendering a detailed and clear picture of the 
dissipative behavior on all time scales for arbitrary strong 
coupling to the leads. 

While previous works have studied the non-equilibrium 
dynamics in the long time limit [Til . l36l . |39| . we have 
provided new information in the intermediate time do- 
main, where exact analytical results are not available. In 



the nonadiabatic limit for strong system-lead coupling 
we have discovered a new non-equilibrium regime with 
a Marcus-like spin relaxation rate. Here, while the non- 
equilibrium dynamics is qualitatively similar to the equi- 
librium dynamics at a finite temperature, the analogy 
is not complete. In particular, a simple linear mapping 
between temperature and bias voltage does not exist, 
in contrast to the electrically damped harmonic oscil- 
lator model [H, HH. The Marcus- like relaxation rate 
exhibits highly nonlinear current-voltage (I-V) charac- 
teristics: The current is practically suppressed at small 
bias voltage, is strongly enhanced at intermediate bias 
(of the order of the energy difference between spin levels 
B), while for large bias linear I-V behavior emerges. 

In the long time limit a non-equilibrium situation gen- 
erates complex scattering phase shifts which are reflected 
in the dynamics through different effects: (i) onset of an 
exponential decay for the spin polarization, (ii) appear- 
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ance of a power law term in the relaxation dynamics, 
with a complex exponent, and (hi) the possible existence 
of an anti-orthogonality-regime. The effects presented in 
this paper are not limited to the specific model utilized 
here, but can be rederived for other systems, e.g. the 
resonant level model of Appendix A, where the polar- 
ization of a spin impurity couples to the resonant level 
occupancy [llj |. 

Going beyond the nonadiabatic limit, we have stud- 
ied multiple spin-flip events with the aid of exact nu- 
merical calculations. Interestingly, we have found that 
the Anderson- Yuval-Hamann treatment of the equilib- 
rium Kondo effect [62j can be extended to the out- 
of-equilibrium regime, but only for weak to interme- 
diate system-bath couplings does the standard pair- 
wise Coulomb gas behavior hold qualitatively for all 
timescales Deviations occur precisely in time inter- 
vals where the Gaussian decay of the correlation function 
Cf(t) is prominent. 

Several future directions are worthy of investigation. 
First, we have restricted ourselves in this work to zero 
temperature. Including the effect of finite temperature is 
straightforward both analytically and numerically [54j . 
In particular, the mapping t — > tanh(7r kbTt) [4lT | trans- 
forms all analytical expressions to those valid at finite 
(but low ksT < Afi) temperatures. Similarly, the nu- 
merical approach of section IV. A may be generalized to 
arbitrary temperatures. The simple model studied here 
can be extended in several important ways, including 
coupling of the quantum subsystem to vibrational de- 
grees of freedom. 

Lastly, we have delineated the precise set of regimes 
where the standard Anderson- Yuval-Hamann Coulomb 
gas behavior is quantitatively accurate. This lays the 
groundwork for future exact numerical studies of the spin 
dynamics for the models discussed here. In particular, 
standard influence functional methodology may be di- 
rectly applied in regimes where pair-wise Coulomb gas 
behavior is exhibited [l3T |. In regimes where deviations 
exist, numerically exact Monte Carlo without the pair- 
wise assumption may be performed. Both of these ap- 
proaches are currently being pursued. 
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APPENDIX A: The Spin-resonant-level model 

We present here a variant of the model system @ 
leading to dynamics analogous to Eqs. (f3"4"| - (j3"51) . The 
model was presented in Ref. [ll| for the analysis of 



the generalized fluctuation-dissipation relation for out-of- 
cquilibrium systems. It describes a spin system coupled 
to a spinless resonant level (creation operator <v ) , which 
is itself coupled to two electron baths n = L,R, 

tl — H s + n B +n SB , 

-B A 
H s = jo- z + —a x , 



H 



if) 

SB 



k , n 



j z dU 



k.n 



(I + °z 



(Al) 



Here B and A are the spin parameters, describing the 
energy gap and the tunneling splitting respectively. J z 
reflects the strength of system-bath interaction, and Vk. n 
is the coupling element of the resonant level to the n- 
th electronic reservoir. I is the identity operator. The 
relation of this model to the generic Hamiltonian ((2J) is 
revealed by diagonalizing H B \ and rewriting Eq. (|A1P 
in terms of the new operators as follows, 



H 



if) 



H 



(!) 



k.n 

J. 



SB 2 
&k,n — ^ ^ T]k,n:k' .n' Ck' 7 n' ] 



k.n.k' ,n' 



d= yVfc,nC fc ,„.(A2) 



with 



n,n' = L,R. The coefficients Vk,n and T)k, n ,k',m & re 



given by [llj 



Vk,- 



V , 

£fe X/fc'.m £ fc — e fc /— iS 
Vk,n v k'n 



5k,k'S n 



e k 



+ id' 



(A3) 



where S goes asymptotically to zero. Next we assume 
that the resonant level- lead coupling is a constant, inde- 
pendent of momentum. The phase shifts, complex num- 
bers in non-equilibrium situations, then become 



Sl = atan 
5r = atan 



A sin 2 (6») 
1 -?:Acos 2 (6>)' 

Acos 2 (6>) 
1 +i\sm 2 (6) ' 



(A4) 



Here A is a dimensionless coupling strength and tan(#) 
determines the asymmetry with respect to coupling to 
the L and R sides 



A = 



2J Z 



r, 



S 9 = atan [y/KJf^j , (A5) 



r n = 2irV 2 p is the hybridization of the resonant level 
with the n-th reservoir and p is the reservoirs [Eq. (|A1[) ] 
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density of states. When the system is symmetric, 
Tr, 8 = 7r/4, and we obtain the following relations 



\S 



5l + 5r = atan(A), 
1 



-ln(l + A 2 ), 



s\ = 



2S L 5 R = 



atan 2 A- iln 2 (l + A 2 ) 

atan 2 A+iln 2 (l + A 2 ) 

4 v ' 



(A6) 



At weak coupling the correlation function can be de- 
rived using the cumulant expansion, as done through Eqs. 



*/(*) 



2- 



ln(l + iDt) + iE s t 

1 — cos(A/xi) 



Si(A/i*) - 



A/it 



^[ 7e + ln(AMt)-Ci(A^)]. (A7) 



with E s = DX/n. By following the derivation which leads 
to Eqs. |Q2J), (20]), $26]) and (J27J) for the present case, 
we obtain the spin-resonant level correlation function at 
strong coupling, 

( {Dt)~P A/it < 1 

C/(t) ~ i (Dt)-^ e - KA " 2 * 2 A/it - 1 - 10 (A8) 

[(Dt)-' 3 (A/it)*> , e- rA ' rf A/it>l, 

with the coefficients 



/3 = atan 2 (A)/7r 2 , 
r = ^ln(l + A 2 ), 
atan 2 (A) 



(A9) 



1 

2^ 



atan 2 A 



iln 2 (l + A 2 ) 



We note that the orthogonality-anti-orthogonality tran- 
sition takes place when the exponential {5\ + S^) changes 
sign, at \ ln(l + A 2 ) = atan(A). 



APPENDIX B: Derivation of the classical Marcus 
rate in the spin-boson model 

In this Appendix we derive the classical Marcus behav- 
ior of a two-level system coupled to an oscillator bath, 
and compare the result to the non-equilibrium Marcus- 
like behavior found in the main text. The classical Mar- 
cus result : 45| emerges in the high temperature limit of 
the asymmetric spin-boson model in the nonadiabatic 
regime [l|. The Hamiltonian is given by 



H = Hc 



H 



H 



(b) 
SB' 



(Bl) 



where the spin system H$ includes a two-level system 
(TLS) with a bare tunneling amplitude A and a level 

splitting B. The reservoir Hg includes a set of indepen- 
dent harmonic oscillators, and the system-bath interac- 
tion H^ B is bilinear in the reservoir coordinates and the 
spin polarization 



Ha 



H 



(b) 



B A 

X, 



(B2) 



Here bj , bj are bosonic creation and annihilation opera- 
tors, respectively. In the nonadiabatic regime the excita- 
tion rate can be calculated within Fermi's Golden Rule 
as [44J 



rf = = 



,±iBt 



C b (t); C b (t) = e 



Mt) = [0- + 2n j )-(l + n j )e 

duj J(w) 



oo 



fljC 



ujf 



4n(w) sin — 



(B3) 



with spectral function J(ui) = tt^j XjS(uj — u)j). Here 
n(u>) 



p u)/k B T 



is the Bose-Einstein distribution 



function with T as the temperature of the bosonic reser- 
voir, Ub is the Boltzmann constant. The case relevant to 
the present paper is the Ohmic spectral density: a con- 
tinuum of bosons with J(uj) — 2-koloj at low frequencies 
and J(u>) — > for u> greater than a cutoff scale u> c . A 
conventional choice is 



J{ui) = 2irauje 



(B4) 



but most of the results do not depend on this choice. It 
is useful to decompose $h into two contributions 



$ fc (t) = + <f> 2 {u c t;Tt). 



(B5) 



Here 



f-cWM sin 2 f 

2 J Q 7T LU 2 e^T-l 



(B6) 



The first term gives the bosonic analogue of the zero tem- 
perature power law dependence ln(l + iDt) found in the 
main text, 



$i(w c t<l) — iXt; X = A 2 /ujj = 2auj c 

j 

$i(w c t>l) — 2aln(iw c f), 



(B7) 
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where A is the solvent reorganization energy. The second 
term gives the analogue of the exponential/Gaussian be- 
havior. For k B T > u) c , the canonical Marcus result is ob- 
tained. In this limit one approximates e u l kBT — 1 — > -^-^ 
leading to 



Marcus 



4fc B T 



duj J{u) sin 2 ^ 



(B8) 



At oj c t -C 1 we approximate sin' 



2 ut 



us 2 t 2 /A and obtain 



•I- 



\ujct < 1) = k B TXt 2 



(B9) 



For uj c t 1 we may set J(w) = 2-Kctoj and get a linear 
behavior, 



Marcus 



(uj c t > 1) = 2irak B Tt. 



(BIO) 



The crossover scale is t* = 1/lv c , with the value 
$f arcMS (i*) PS 2ak B T/uj c . Thus, for any a, at large 
enough k B T/uj c , <& 2 becomes large enough that a Gaus- 
sian relaxation results. 

The analogy between the high temperature Marcus 
behavior and the results we have found in the non- 
equilibrium fermionic model is not complete, since in the 
latter case we typically assume that the electron bands 
are wide relative to the potential bias. We therefore con- 
sider next the analogous equilibrium limit of k B T < to c . 
In this case the short time limit (sin(wi/2) — * Lot/ 2 ) of 
Eq. (HU) obeys 



§l{uj c t < 1) = 2at 2 



duj- 



j/u 



o/k B T _ l 



at 2 k 2 B T 2 ir 2 



(Bll) 



while the long time limit reduces to (|B10[) . We thus 
obtain a short time t 2 and a long time i-linear behav- 
ior. The crossover occurs at t* ~ 1/fegT and $^(t*) = 
7r 2 a/3 ~ 3a. Thus, for a much smaller than 1, the re- 
laxation integrals are dominated by the long time region 
where $ ~ t, leading to an exponential relaxation. For 
a > 0.5 however, the power law prefactor ensures that 
the integral is dominated by short times, of order lo~ x , 
so that the frequency dependence is significant only on 
the scale of the cutoff scale u> c . 

We compare next this behavior to the non-equilibrium 
Marcus-like behavior found in the main text. The math- 
ematical essence of the non-equilibrium result is that the 
time decay function behaves as <E>/ ~ n(A^t) 2 at short 
times, and as $/ ~ TA/it at long times; the crossover be- 
tween these two regimes occurs at t* ~ T/kA/j, and the 
value at t* is <S>f(t*) ~ T 2 /k. In the weak coupling limit 
r ~ k, t* ~ 1/A/x and $/(<*) ~ T 2 / n ~ v % < 1, so the t 2 
behavior is not important for the relaxation rates. How- 
ever, in the strong coupling limit, T/n » 1, t* 1/A/i 
and $/(t*) 3> 1, so that by the time t reaches t* the evo- 
lution function has become negligibly small. In this cir- 
cumstance the t 2 behavior controls the relaxation (for all 



relevant energy differences), leading to the Gaussian be- 
havior discussed in the text. In the non-equilibrium case, 
the key parameter is therefore r 2 /K, and as this becomes 
larger than unity, Gaussian behavior results. In contrast, 
in the classical (high temperature) Marcus limit, the role 
of T 2 /k is replaced by 2ak B T '/u c , while for k B T < ui c , 
the role is played by a. If this is small, one has expo- 
nential relaxation, while if this is larger than 1/2, the 
kinematics are different and the frequency dependence 
is controlled by the bandwidth scale lu c . Thus the non- 
equilibrium Marcus-like rate found here is really a new 
phenomenon. 

We proceed and calculate the classical Marcus rate in 
the high temperature limit. We substitute Eq. (1B9|) and 
the short time limit of Eq. (|B7j) into Eq. (|B3|) . perform 
the Fourier transform, and recover the Marcus relation 
for the nonadiabatic rate 



D -(TB+\) 2 /(4\k B T) 



\k B T 



(B12) 



The temperature dependence of the rate constant shows 
an activated regime for |A ±-B| ^ 0, while in the absence 
of the barrier, — B = A, the rate decreases with T. The 
excitation rate and the emission rate rjj" are related 
one to another through an activation factor as 

i7/r+ = e -w. 



(B13) 



For completeness, we include here other results of the 
spin-boson model in the nonadiabatic limit: The Golden 
Rule rate to lowest order in k B T/u c and B/lj c yields 
(B>O)0 



A 2 / LO c 

4^~ c \2irk B T 
\f(a-iB/2Trk B T)\ 2 _ B 



/2k B T 



T(2a) 



(B14) 



Here T(x) is the complete Gamma function. For weak 
damping, a -C 1, this expression reduces to 



2a naA 2 



eff 



2nk B T 

A eff J [( 2n 7r/.-i,. T Y- + B- 
B 



where A e // is an effective tunneling element {a < 1), 

A e ff = A[f(l - 2a)cos(7ra)]^T (A/uj c ) a ^ 1 ~ a '> . 

(B16) 

At zero temperature we can calculate the rate exactly for 
an arbitrary cutoff frequency (B > 0) 



r+(T = o) 



ttA 2 



2T(2a)uj c \uc 



B 



-B/u 



rr(T = o) = o. 



(B17) 



As expected, at T = there are transitions only from the 
upper level to the lower state. 
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APPENDIX C: Bosonization of the non-equilibrium 
Fermi edge Hamiltonian 

We briefly present here some of the relations be- 
tween bosonic and fermionic operators, and transform 
our fermionic system-bath Hamiltonian into its bosonic 
analog via bosonization p6l |. and discuss the issues in- 
volved in bosonizing the non-equilibrium version of the 
model. For simplicity, instead of the spin-boson Hamilto- 
nian Eqs. (U)-©, we discuss here the x-ray edge Hamil- 
tonian (d2), describing the interaction of a localized core 
hole with two (possibly out-of-equilibrium) metal leads 



k,n—L,R 



k,k f ,n — L,R 



+ Y \v L ,Ra\ L a k ^ R + V R xa\ R a k ^ L S d. (CI) 

k,k' 

The first term includes two isolated Fermi baths. The 
second and third terms describe intra-bath processes 
(diagonal coupling), and inter-bath interactions, respec- 
tively, d (d^ ) are annihilation (creation) operators of the 
core hole. We consider bands of width Do and a constant 
density of states p. 

To solve the equilibrium problem one proceeds as fol- 
lows. First, one defines new fermion operators 



£*fc,+ = cos-a k ,R + sm-ak,L 



(C2) 



u k - = -sin-a k ,R + cos -a k ,L (C3) 



with 



+ 9 V RR - V LL 
tan - = — ■ 

2 2^V L rVrl 

in terms of which Eq. (|C1[) becomes 



(C4) 



k,n— + , — k,k' ,n— +,— 

(C5) 



with V± = v ^+ v ™ ± V {Vll+ ^ rr)2 + VlrVr L . Cru- 
cially, in equilibrium the new fermion variables obey the 
usual Fermi statistics 



S k , k 'S n ^ n 'f(e k /k B T), 



(C6) 



with /(ek/ksT) as the Fermi-Dirac distribution function. 
For this reason, each channel can be bosonized, leading to 
the standard result of Schotte and Schotte for the Fermi 
edge singularity problem 25]. In a non-equilibrium situa- 
tion, while the change of basis can be made, the different 
distribution functions for the left and right leads mean 
that Eq. ([C6j does not hold, preventing bosonization in 
the transformed basis. One may attempt to proceed by 
defining the density operator 



»(?) = Y a k+q,n a k,n: (jl = L, R) , 



(C7) 



which creates particle-hole excitations in the n lead with 
momentum q. We use it and define bosonic creator and 
annihilator operators that obey the bosonic commutation 
relation 



bin = \ljZPn(l) ( 'i • (,) - 

2T 



b q ,n = \l —p n (-q) (q > 0). 



(C8) 



The distance L is related to the density of states through 
p = L/2i:vf with vf as the velocity at the Fermi energy, 
taken to be the same for both reservoirs. The fermion 
field operators 



= 4f Yl afe .«> *n = -7f J2 a L 



(C9) 



can be expressed in terms of the boson operators as [2 



: exp 



/27T 

^VTq e ~ 
q 



aq/2 



{bq.n b Qn ) 



(CIO) 

Here F n (F£) are the Klein factors that lower (raise) the 
total fermion number in the n reservoir by one. The 
chemical potential difference is therefore concealed inside 
these factors, a is an arbitrary cutoff that regularizes 
the theory and mimics a finite bandwidth. Using these 
expressions, the fermionic Hamiltonian (|C1[) translates 
into a bosonic expression as follows, 



H = H B + (A Q + A b )dU. 
H includes the isolated reservoir term 
H B = Y v Fqbl tn b q , n , 

q,n 



(Cll) 



(C12) 



and diagonal (A a ) and nondiagonal (Ah) contributions 



n.n \ 



q . n 



Ai, 



(V l ,rU l U r + Vr, l UrUI). (C13) 



Here 



Un 



j? ne -E,A,(6+, n -6,, n ) j ( n = L,R); 
VTq e ' 



*q/2 



(C14) 



All the prefactors are absorbed into the coefficient 
Vl,r = ^V l ,r. Eqs. (fCTT]) - (fC13|) reveal that for a 
non-equilibrium system bosonization yields a nonlinear 
Hamiltonian with a highly nontrivial form. Compared 
with the solution in the fermionic picture (Appendix D), 
bosonizing the Hamiltonian does not simplify the calcu- 
lation, as it does in equilibrium. It should be noted, how- 
ever, that the use of Eqs. (|C11|) . (|C13j) and (|C14|) yield 
second cumulant expressions identical to those derived in 
Appendix D. 
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APPENDIX D: Derivation of the second cumulant 
expression 

We derive here the details of the weak coupling corre- 
lation function Eq. (|2"2")l . The second cumulant is given 
by 

K 2 {t) = -\l dh [ dt 2 (TF(t 1 )F(t 2 )) c , (Dl) 
z Jo Jo 

where F = J2k,k>,n,n' V n,n'%,n a k',n'> (...) c denotes a cu- 
mulant average and T denotes time ordering. For simplic- 
ity, we disregard diagonal interactions, Vl,l — Vr.r = 0. 
The integrand is calculated using Wick's theorem [44| to 
yield 



C 2 (r) = (TF(t)F(0)) c 







(pVl,r) 2 


/ de 1 




J-Da 


ff'R 




/ dea 


/ de\e~ 


J-D a 





Do 

de 2 e l{£1 - e2)T 

mi 

i(ei-£2)T 



(D2) 



Here 2Z?o is the bandwidth, p is the energy independent 
density of states, and pl (pr) is the chemical potential 
at the L (R) lead. This expression assumes zero temper- 
ature. Assuming wide bands, pk < D$, the integrals in 
(ID2[) can be trivially performed producing 



C a (r) = -2( P V L , R ) 2 

D t>1 



1 



-iD r \ 2 



2(pVi 



L.R) 



D 2 



(1 + 1D t) : 



cos(A/rr) 



■cos(A/zt), (D3) 



with the voltage difference A/i = pl^P-r- In equilibrium, 
the correlation function is therefore given by 



c 2 q (t) 



'{i + iD ty 



(D4) 



We substitute this expression into Eq. (|D1|) . and obtain 
the second cumulant approximation for equilibrium situ- 
ations 



Kl\t) = -{pV L , R ) 2 \n{l+Dtf). 



(D5) 



This is the standard result for Tomonaga's model 25]. 
In non-equilibrium situations, A/i ^ 0, the correlation 
function includes an oscillatory function, (|D3|) , which can 
be decomposed into its equilibrium and non-equilibrium 
contributions as follows: 



2(pV LiR ) 2 D 2 



(1 + iD ty 



[cos(A/ii) - 1] . (D6) 



The equilibrium term yields Eq. (|D5[) . We proceed with 
the non-equilibrium part. For A/i <C Dq, D$t > 1 

Jo Jo 

\2 , 



-( P v L , R yAfi / dh / rft 2 

Jo Jo 



ipVi 



L.R) 



dh I dt 2 — 



{h-t 2 y 

sin(A/i(ii - 1 2 )) 

h-t 2 

cos(A/i(ti - t 2 )) - 1 



ti-t 2 



(D7) 



Exact integration leads to 



-2{pV L ,R.) 2 [AptSi(Apt) - (1 - cos(A//f))] 
2(^ L , fl ) 2 [ 7e + ln(A/ii) - Ci(A M t)] . (D8) 



The sine and cosine integrals are defined as Si(x) = 
Jo* sJ ^dt, Ci(x) = 7e + Hx) + /* and 7e = 

0.5772 is the Euler-Mascheroni constant. In deriving 
(|D8[) we have used the following identities: / Si(a;)dx = 
cos(x) + xSi(x); /* sin[{ ^ m dt = Si[(x - a)j3] - Si(a/3). 
The sum of Eqs. (|D5[) and (|D8|) is our expression for the 
second cumulant (|22[) with = tt/oVi,^. After exponen- 
tiating, the first term provides an exponential relaxation 
at long times, while the second term yields a power law 
contribution. 



APPENDIX E: Approximate Analytical Evaluation 
of Rate Constants 

In this Appendix we present a derivation of the rate 
constants in the important limiting cases. We begin from 
the basic expression (|33|) 



rf(o;) =5R 



dt 



(iDt) 



0effW 



a iuit — A(t) 



(El) 



Here D is an energy scale of the order of the bandwidth, 
and P e f f is an effective exponent which changes from the 
equilibrium power (3 to the non-equilibrium power (3 neq . 
The results of section IIVI imply that 



A(t) 



(E2) 



with <j) defined such that 4>{x — > 0) — > x 2 and <fi(x —> 
oo) — * x. This implies that the coupling constant g and 
characteristic time t* are 



K 



nAp 



(E3) 
(E4) 



In the weak coupling limit, (y <C 1), T ~ k ~ v 2 so 
g <C 1 and i* ~ (A/i)" 1 , while in the strong coupling 
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limit {y ~ 1), Y — > oo while k saturates, so j > 1 and 
t* 3> (A/i) _1 . We define a dimensionless time coordinate 
u = t/t* and frequency oj* — ojt* in terms of which the 
dimensionless relaxation rate DT^ becomes 



r- 

DT+(oj) =5R / 
Jo 



(iu) 



(E5) 

The analysis of the integral in Eq. (|E5|) requires some 
care because the scales which dominate the integral may 
not be the scales which dominate the real part of the 
integral. To isolate the contributions to the real part, 
we deform to contour into the complex plane. Writing 
u = x + iy we deform the integration contour into two 
parts, one running along the imaginary axis (x = 0) to 
the point y = y* at which ioj* — gd<p/dy = and another 
running parallel to the real axis along the contour u = 
x + iy* . Thus we have 



DT+{oj) = h+I 2 , 



with 



Ji = 9f 



(E6) 



(E7) 



— e -\u y 



5R 



(22; _ y*\0eff(x+iy*) 
iuj*x-gtf>(x+iy*) (E8) 



Here 3? refers to the real part of the integral and 3 to the 
imaginary part. We analyze these equations first in weak 
coupling g < 1. Let us begin with I-y. Inspection of the 
second cumulant formula shows that y* < if oj < 0. In 
this case the integral does not have any imaginary part 
so for u> < 0, 7i = 0. Next consider small positive u, 
where we may approximate <fi(y) — y 2 implying 



V" = 



2g 2TAfi 



(E9) 



Thus for oj < rA/j we may use the equilibrium exponent 
and approximate <f> — y 2 obtaining 



h = Q{oj) sm{n f3)(Dt*) 



^ dy 

7 ( 



-v y+gy 



(E10) 

At the endpoint ^- the argument of the exponential is 
minimized; the minimum value is 



4g 4fcAp 2 



Substituting the maximum value u> = TA/i and noting 
that in weak coupling r ~ n -C 1 we see that the argu- 
ment of the exponential is negligible over the entire range 
oj < TA/i and we get 



h = e(w) 



sin(7r/3) / Doj 



1-/3 \2TAn 2 



i-0 



(Ell) 



For oj > TA/i, y* > 1, and we must consider the form of 
<f> for large imaginary argument. Inspection of the second 
cumulant formula shows that 



(j>(iy) ~ cosh(y) ~ e v , 



(E12) 



implying y* ~ \xi(u)*/g). In this case the value of the 
argument of the exponential at the upper limit of inte- 
gration is -oj* 1ii(oj* / g) ~ — (w/A/i) ln[w/(rA/x)]. Thus 
for w > j^yp the upper limit of the integral may be set 

to infinity, and for oj > A/x the integral is dominated by 
y < 1 yielding 



Drt(w) - Q(-l)- /5 (Dr) 1 -' 3 / ^e^* z (E13) 

Jo zl 

D ^-0 



sin(7r/3) 



r(l-/3). (E14) 



Here f (x) is the complete Gamma function. Eq. (|E14|) is 
simply the usual equilibrium result. Therefore, in weak 
coupling, Ii is given by Eq. (|E11[) for oj < A/i/lnT -1 
and by Eq. (|E14|) for cj > A/i, with a rather broad 
crossover regime. 

We next turn to 1% which is non-zero for both signs 
of oj. The frequency regimes are as for I\. For oj < 
A/x/lnT -1 the integral is dominated by large x where 
4> — x and the prefactor e~' u v ' is negligible, so that we 
find 



h « {D) 1 -^T{l-p neq )St 



-i7T/3 Tle ,/2 



TA/z) 1 



-/3„ 



(E15) 

Note that we have replaced (3 by the long time non- 
equilibrium value Pneq- In the weak coupling limit, 



(3 n , 
behavior 



eq 



< 1. Setting /3 



neq 



yields a Lorentzian 



DTAfj, 

! + r 2 A ) u 2 



(Ei6) 



Note that for w > I2 only becomes smaller than I\ for 
w ~ rA^/sin/37r ~ A[i. As w becomes of the order of 
A/x/lnT -1 the prefactor begins to be important and I2 

decays proportional to e^ 3 ^ as found by Mitra et 

al [3. 

To summarize, for weak coupling we find a rate which 
for small frequencies oj < A/i/lnT -1 is approximately 
Lorentzian, with decay constant TA/i. On the emission 
(positive frequency) side, the Lorentzian decay is over- 
come by the contribution of I\ and eventually crosses 
over the equilibrium rate, Eq. (|E14[) . while on the ab- 
sorption (negative frequency) side the rate crosses over 
to into the e - ^ ln (3jr) A/x relaxation. 

We now take up the strong coupling (<? > 1) limit. 
For oj < rAp Eq. (|E10|) still applies, but now at the 
endpoint of the integration region the argument of the 
exponential can be large. For oj* < yjg (i.e. oj < ^/kA/i), 
the variation of the exponent is not important, and we 
get 



^) =0 (^ Sin(7r/?) 



1-/3 \2nAfi 



Doj 



1-/3 



\oj\ < y/HAfi. 

(E17) 
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In contrast, in the opposite high frequency limit we can 
set the upper limit of the integration to infinity and drop 
the y 2 term. This yields 



l[ B) =9(o;)sm(7r/3) ( — 



1-/3 



f(l-/3); 



uj > ^/kA^i, 
(E18) 



which is again the equilibrium (A[i = 0) result. This 
continues to apply even for uj > TA/i. We next turn to 
1%. For uj < TAfi we again approximate (j) = (x + iy*) 2 
and find 



7 2 _ e'^^iDt*) 1 - 13 ^. 



dx 



-gx 



(E19) 

The important x are of order Jg so that for uj < y/nAfi 
we may neglect the u> in the denominator and get 



r (A) 



7T/3 

cos — e 
2 



r [(1 



/3)/2], 



In the opposite limit |w| > ^fnAfi we neglect x in the 
denominator and find 



T (B) 



^ cos{-k(3) 



D 



2kAh 2 
2^/k.A^ \ Duj 



\uj\ > v^kA/z. 

(E21) 



For the emission rate uj > 0, the Marcus rate given by Eq 
IE21l goes over to the equilibrium power law behavior, Eq 



IE181 when 1^ ' becomes smaller than /} ' , when happens 
for uj slightly larger than 2y / KA/i. 



Finally, if uj* > g (uj > TA/z) then the approximation 
< v^A/x.0 = y 2 does not apply and the rate goes over the e _wlncJ 
form discussed above. However, by this time the rate is 
(E20)so small that this behavior is not relevant. 
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